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Abstract 

The qualitative aspects of the phase diagram of the Ising model on the cubic lattice, with ferro- 
magnetic first neighbors (Ji) and antiferromagnetic second neighbor couplings (J2) are analyzed 
in the plane temperature versus a, where a = J2/J1 is a frustrated parameter. We used the orig- 
inal Wang-Landau and the standard Metropolis algorithm to compare past results of this model 
obtained by the effective-field theory for the cubic lattice. Although the nature of some critical 
points, chosen at relevant values of a, show that the phase diagram is, in general, qualitatively 
correct, our Monte Carlo results suggest that the reentrance form of the frontier that separates 
the ferromagnetic and the colinear order is an artifact of the effective-field theory, which might 
disappear by improving these approach. 
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I. INTRODUCTION 



Some magnetic compounds like Eu x Sri_ x S [l, 2| and Fe x Zn\- x Fi {3] present more than 
one low-temperature magnetic orderings, depending on its parameters, like the strength of 
the interactions and the concentration of magnetic ions x. They are well- described by mod- 
els which consider competitions of nearest-neighbor and next-nearest-neighbor interactions. 
The simplest model which may describe such compounds is represented by the following 
spin hamiltonian: 



H 



(1) 



where <Ji = ±1, and % = 1, . . .,N, where N is the total number of spins. The first sum 
contains all pairs of ferromagnetic nearest-neighbor couplings (Ji > 0), and the second one 
is for all the next-nearest-neighbor antiferromagnetic interactions (J2 > 0). The model with 
J2 < 0, is well understood and establishes the Ising second-order universality class J4]. Nev- 
ertheless, this model has attracted a lot of interest in the past, especially when implemented 
in the square lattice 5M29] . For this case, the magnetic order at zero temperature depends 
on the value of the frustrated parameter a = Jil\J\\- For a < 1/2, the order is ferromag- 
netic (F, Ji > 0) or antiferromagnetic (AF, J\ < 0), and for a > 1/2, we have the collinear 
order, also called superantiferromagnetic order (SAF). So, for 1/2 < a < 1, there has been 
controversial results about the nature of the order-disorder transition at finite temperatures. 
Recently, Kalz and Honecker [30J have concluded that Monte Carlo (MC) data, obtained 
with lage sizes (L = 1000, 2000), yield a clear picture only for 1/2 < a < 1, where a first- 
order phase transition scenario is established by the double peaked structure of the energy 
histograms, and for a > 1 there must be continuous phase transitions. 

In the present work we study the model in Eq. ([1]), implemented on the cubic lattice with 
periodic boundary conditions, and J\ > 0, J2 > 0. 



'his model has already been treated 
within an effective-field theory by dos Anjos et al. [31]. So, in Fig. ([1]), we show the phase 
diagram of this model in the plane fc^T/Ji — a, recalculated by the present authors using 
effective-field theory with a cluster of one central spin (EFT-1) as done in reference [31] . At 
zero temperature, it can be exaclty determined two type of orderings separated by a = 1/4. 
For a < 1/4, the ferromagnetic order appears, whereas for a > 1/4 a SAF order is set. 
At finite temperatures those phases are separated by a first-order frontier, which presents a 
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reentrant form, as shown in the inset of Fig. ([I]). This frontier finishes at a critical end point 
(CE), where two order- disorder frontiers are also ending. The first one is of second-order 
type and separates the F-Paramagnetic (P) transition, for lower values of a, and the second 
one is of first-order type and it separates the SAF-P transition, for higer values of a. 

In that work, the authors used a decoupling procedure, which ignores all high-order corre- 
lations so as to approach the unmanageable expressions of all boundary spin-spin correlation 
functions. Although this analytical treatment improves the mean-field approach, which is 
insensitive to frustration, accuracy and qualitative aspects can be lost in determining the 
critical temperatures and the nature of the phase transitions. In order to verify the quali- 
tative aspects of these phase diagrams, we need to use powerful Monte Carlo techniques to 
construct the canonical probability distribution function (CPDF) (P(E,T) ~ exp(—/3E)), 
for given temperatures, and finite sizes of the cubic lattice. Accordingly, at a critical tem- 
perature the CPDF will show a double-peaked form for a first-order phase transition, and 
a single-peaked form, for a second-order one. So, the original Wang-Landau sampling al- 
gorithm (WLS) is a suitable MC method to get the CPDF from density of states g(E) 



32 



33|. 



One of the advantage of this method is that we directly construct the density of states 
g(E,T), through which the canonical partition function is achieved, so all the thermody- 
namic variables can be plotted as a temperature's function (free energy, heat capacity, etc). 
Furthermore, at low temperatures, the Metropolis algorithm will get trapped in states of 
energy local minima at low temperatures 34j, especially in frustated models. For instance, 
conventional simulations in the canonical ensemble would not be efficient in the region close 
to a = 1/4 (see Fig. (pQ)), where the system is in a highly frustrated zone. Nevertheless, the 



original WLS is not without accuracy problems j35|, however, it does not affect qualitative 
results. Another problem appears when larger lattice sizes are needed. In this case, we 
require to divide the relevant energy range into fixed windows, then we have to join them 
after convergence is achieved. Consequently, the resulting density of states and associated 
thermodynamic functions are shown to suffer from boundary effects. This undesirable effect 
becomes more conspicuous for the obtention of g(E,M), which is useful to calculate the 
canonical probability distribution function including the order parameter P(E, M,T). In 
this case, it is necessary to perform a two-dimensional random walk in a relevant (E, M) 
space. In most cases, this relevant (E, M) space needs to be divided by surfaces, but after 
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matching them the resulting P(E, M, T) will have small discontinuities. 

The general source of these difficulties seems to be due to the difficulty in matching 
surfaces at the boundaries rather than curves as in one- dimensional random walks ^M^. To 
overcome this problem Cunha-Netto et al. proposed the WLS with adaptative windows {37) ]. 
where instead of defining fixed energy windows, the boundary positions depend on the set 
of energy values on which the histogram is flat at a given stage of the simulation. So, errors 
that may arise near the border of a given window are corrected in subsequent stages, in 
which the border positions are shifted. Nevertheless, it improves the quality of the results 
of the WLS with fixed windows at the expense of computational cost. The WLS algorithm 
with adaptative windows considerably increases the computational time with the size of the 
system, and seems not to be able to be parallelized. Therefore, in this paper we use the 
multi-range original WLS algorithm with fixed windows, which does not affect qualitative 
results as will be shown. 

II. METHODOLOGY 

We used the original Wang-Landau algorithm so as to get the corresponding logarithm 
of the density of states (log g(E)) for the model defined in Eq. (pQ). Consequently, we 
can calculate the mean energy E and specific heat curves C, and the canonical probability 
distribution functions P(E) therefrom. These are our least necessary tools to do a qualitative 
analysis of the criticality of the system. In order to get log g(E), the minimum E min and 
maximum E max energies of the system are needed, for a given value of a and L. Then we 
also need to number every discrete energy value Ej between them to define an integer array 
H(Ej) and a real array g{Ej) as useful histograms for the algorithm. Initially, the g(E) is 
unknown, so all bins in the array are set to unity. Since the typically range of g(E) is of 
high orders of magnitude, it is common to store log g(E). In addition, a visits histogram 
H(E) is maintained. 

Initially, all bins have zero visits for both log g(E) and H(E). The bins are then filled over 
the course of a MC simulation, in which moves (spin flips) are accepted if p < min < 1, ^r^k \ , 



where p is a uniform random number in the range [0,1], E and E 1 are the energies of the 
current and the proposed move, respectively. After the move is accepted or rejected, the 
histogram H(E) is incremented by one and the density of states histogram g(E) is multiplied 
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by a constant factor /, such that g(E) =— > g(E) x /, where the initial choice is / = e ~ 2.72. 
An accurate estimate of g(E) is reached if the histogram H(E) becomes flat. 

At this step the histogram H(E) is set to zero and the modification factor / is reduced 
such that f i+ i — > tfjl. This process is repeated until fa be close to 1, so we repeat it 
until % = 14, using n = 4 to accelerate the process. However, the repetition of the above 
simulation suffers from the shortcoming that very large entries need to be stored in g(E). As 
mentioned before, in order to avoid this problem, the quantity log g(E) — > log g(E) + In/, 
is evaluated. The modification factor is then now updated as log(/j+i) — > (1/n) ln(/j). 

The adopted flatness criterion was H(Ej) > 0.8 (H(E)), V j. However, for the present 
model, it is difficult to satisfy it around a = 0.25, on account of frustration. So, for a 
given fi, we stop the process after a maximum number of Monte Carlo moves (M max ). On 
the other hand, it is important to mention that it is not necessary to use the entire energy 
interval [Emm ■ ■ ■ E max ] of the system to get the relevant information of the criticality. Thus, 
we need just to obtain the density of states for the relevant energy subspace [E\, E 2 ] in order 
to calculate the thermodynamic quantities throughout the temperature range of our interest. 
For our model, and for a given lattice size L, the number of energy bins are considerably 
increased for some values of a, this is why we have to apply a multi-range Wang-Landau 
algorithm with fixed windows even for the relevant energy subspace [E\, E 2 }. Otherwise, the 
flatness criterion will never be satisfied. 

III. RESULTS AND DISCUSSION 

We study the model defined in Eq. ([[]) by performing WLS and Metropolis algorithm for 
< a < 1, for the cubic lattice with N = L x L x L sites. We choose L = 16, because for 
L > 16 the number of energy bins are considerably increased for certain values of a. So, too 
many windows would be necessary to apply the multi-range WLS algorithm, which would 
also increase the computational cost. In Fig. (2a), it is shown the logarithm of the density 
of states log g(E), for a = 1.0, for the entire energy space. This is an assymetric function 
in E, in contrast to that of the simplest spin-1/2 Ising Model. 

Figure (2b) shows the corresponding mean energy and the specific heat versus tempera- 
ture obtained therefrom. In Fig. ([3]), we show the CPDF for three different temperatures 
for a = 1. A double-peaked structure appears at the estimated pseudo-critical temperature 
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showing a first-order phase transition. In Fig. (j4J) it is shown the energy range used to 
perform the WLS process for a = 0.5 and for a given step of the algorithm (i — 3). There 
we may see the results for the four overlapped windows in which the selected energy range 
was divided. Then we meet the curves of the four windows into one curve to get the loga- 
rithm of the density of states plus a constant. Accordingly, in Fig. (j5J) we show the relevant 
results for a = 0.5. In figure (5a) we see the specific heat curve in which both, WLS and 
traditional Metropolis simulations are in agreement. Accordingly, it is shown a first-order 
transition as presented in figure (5b). For a = 0.25, the F and SAF orders coexist at T = 0. 
So, WLS results at finite temperatures present a first-order phase transition at the specific 
heat peak as shown in figures (6a) and (6b). It is necessary to point out that for this value 
of a, frustration gets the extreme, so traditional Metropolis simulations are not suitable to 
equilibrate the system, specially at low temperatures. This is why by the knowledge of the 
function of the density of states, obtained by WLS, one may overcome the limitations of the 
traditional Monte Carlo technique. 

On the other hand, in order to verify whether a reentrant behavior occurs for the 
F-SAF frontier, as shown in Figure 1, we explored a closer value of a around 0.25, such 
as a = 0.24. For this case, a relevant energy subspace requires many energy bins even for 
single temperature calculations. Accordingly, we had to construct the specific heat curve 
by sections within the WLS method. In figure (7a) we present these sections, which fit 
well with Traditional Metropolis simulations. At the specific heat peak, the corresponding 
CPDF shows a second-order phase transition in figure (7b), because a single-peaked 
structure appears. Furthermore, we explored the CPDF along the entire relevant range 
of temperatures and no double peak CPDF appeared for a = 0.24, for this considered 
lattice size L = 16. So, this result strongly suggests that there is no reentrance for the 
F-SAF frontier. If that reentrance exist, it must be inside the interval 0.24 < a < 0.25. 
Consequently, we might infere that the reentrance shown in Fig. (CQ) be an artifact of 
EFT-1 approach. On the other hand, we may approximately locate empirically the end of 
the F-SAF frontier by observing the behavior of the orders parameters as functions of a, as 
shown in Fig. (jSJ). So, by the aid of this figure we estimate the location of the CE point 
around k^Tj J\ = 1.0, at a = 0.25. 
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Finally, Fig. (Q shows the phase diagram of the present model in the plane ksT / J\ — a, 
for L = 16. For < a < 0.25 there is a second-order F-P critical frontier. On the other 
hand, for 0.25 < a < 1 there is a first-order SAF-P critical frontier. We affirm that the 
F-SAF frontier seems to be vertical. The three critical curves meet at a critical end point 
which must be around the end of arrow shown in this figure. Therefore, Fig. shows a 
phase diagram which qualitatively agrees with that of the Fig. ([I]) obtained by EFT-1, with 
the exception of the reentrance appeared in it. 

IV. CONCLUSIONS 

The phase diagram of the Ising model in the presence of nearest- and next-nearest- 
neighbor interactions on a simple cubic lattice with N = L x L x L sites, was studied by 
performing MC simulations considering the original WLS and the traditional Metropolis 
algorithm, for L = 16. The transition from the ferromagnetic phase to the disordered 
paramagnetic phase is of second-order type. On the other hand, a first-order transition 
frontier is suggested from the SAF phase to the P one, as well as from the SAF phase to 
the F one. The reentrance appeared in the F-SAF critical frontier obtained by an EFT-1 
seems not to exist for the present model, at least out of 0.24 < a < 0.25. It suggests that 
this reentrance is a consequence of the nature of that approach. However, MC results give 
qualitatively the same phase diagrams as obtained by effective-field calculations. 

ACKNOWLEDGEMENT This work was partially supported by CNPq (Edital Uni- 
versal), CAPES, FAPERJ and FAPEAM (Programa Primeiros Projetos - PPP) (Brazilian 
Research Agencies). 



[1] H. Maletta, in Excitations in Disordered Systems, edited by M. F. Thorpe (Plenum Press, 
New York, 1982). 

[2] H. Maletta, in Heidelberg Colloquium on Spin Glasses, edited by J. L. van Hemmen and I. 
Morgen-stern, Lecture Notes in Physics Vol. 192 ( Springer- Verlag, Berlin, 1983). 



7 



[3] D. P. Belanger, in Spin Glasses and Random Fields, edited by A. P. Young (World Scientific, 
Sin-gapore, 1998). 

[4] R.J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic, London, 1982). 

[5] S. Katsura and S. Fujimori, J. Phys. C 7 (1974) 2506-2520. 

[6] M. N. Barber, J. Phys. A: Math. Gen. 12 (1979) 679-688. 

[7] F. Y. Wu, Phys. Rev. B 4 (1971) 2312-2314. 

[8] J. Oitmaa, J. Phys. A: Math. Gen. 14(1981) 1159-1168. 

[9] R. H. Swendsen and S. Krinsky, Phys. Rev. Lett. 43 (1979) 177-180. 
[10] D. P. Landau, Phys. Rev. B 21 (1980) 1285-1297. 
[11] K. Binder and D. P. Landau, Phys. Rev. B 21 (1980) 1941-1962. 
[12] D. P. Landau and K. Binder, Phys. Rev. B 31 (1985) 5946-5953. 

[13] J. L. Moran-Lopez, F. Aguilera-Granja, and J. M. Sanchez, Phys. Rev. B 48 (1993) 3519-3522. 
[14] J. L. Moran-Lopez, F. Aguilera-Granja, and J. M. Sanchez, J. Phys. Cond. Matter 6 (1994) 
9759-9772. 

[15] C. Buzano and M. Pretti, Phys. Rev. B 56 (1997) 636-644. 

[16] K. Tanaka, T. Horiguchi, and T. Morita, Phys. Lett. A 165 (1992) 266-270. 

[17] J. A. Plascak, Physica A 183 (1992) 563-573. 

[18] P. M. Oliveira, C. Tsallis, and G. Schwachheim, Phys. Rev. B 29 (1984) 2755-2760. 

[19] H. W. J. Blote, A. Compagner, and A. Hoogland, Physica A 141 (1987) 375-402. 

[20] M. P. Nightingale, Phys. Lett. A 59 (1977) 486-488. 

[21] H. W. Blote and M. P. Nightingale, Physica A 134 (1985) 274-282. 

[22] P. A. Slotte, J. Phys. C 16 (1983) 2935-2951. 

[23] H. J. W. Zandvliet, Europhys. Lett. 73 (2006) 747-751. 

[24] F. Aguilera-Granja and J. L. Moran-Lopez, J. Phys. Cond. Matter 5 (1993) A195-A196. 
[25] A. Malakis, P. Kalozoumis, and N. Tyraskis, Eur. Phys. J. B 50 (2006) 63-67. 
[26] J. L. Monroe and S.-Y. Kim, Phys. Rev. E 76 (2007) 021123-1021123-5. 
[27] A. Kalz, A. Honecker, S. Fuchs, and T. Pruschke, Eur. Phys. J. B 65 (2008) 533-537. 
[28] Rosana A. dos Anjos, J. Roberto Viana, and J. Ricardo de Sousa, Phys. Lett. A 372 (2008) 
1180-1184. 

[29] A. OHare, F. V. Kusmartsev, and K. I. Kugel, Phys. Rev. B 79, 014439-8(2009). 
[30] A. Kalz and A. Honecker, Phys. Rev. B 84, 174407(2011). 

8 



[31] Rosana A. dos Anjos, J. Roberto Viana, J. Ricardo de Sousa and J. A. Plascak, Phys. Rev. 

E 76, 022103(2007). 
[32] F. Wang and D.P. Landau, Phys. Rev. Lett. 86, 2050(2001). 
[33] F. Wang and D. P. Landau, Phys. Rev. E 64, 056101(2001). 
[34] Y. Okamoto, AIP Conf. Proc. 690, 248(2003). 

[35] A. A. Caparica, A.G. Cunha Netto, Phys. Rev. E 85, 046702(2012). 
[36] Shan-Ho Tsai, F. Wang and D.P. Landau, Braz. J. Phys. 38, 6(2008). 

[37] A. G. Cunha-Netto, A. A. Caparica, Shan-Ho Tsai, R. Dickman, and D.P. Landau, Phys. 
Rev. E 78, 055701(2008). 



9 




FIG. 1: Phase Diagram of the model described by the hamiltonian in Eq. (fTJ). The vertical axis 
represents the reduced temperarute, and the horizontal one represents the parameter of frustration 
(a = Jij J\). Dashed and continuous lines represent first and second order critical frontiers, 
respectively. The inset shows better the reentrant form of the frontier separating the F and SAF 
orders. CE is the critical end point shown. 
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FIG. 2: (a) Logarithm of the density of states for a = 1 and size L = 16, obtained by the original 
WLS algorithm applied for the whole energy range, without dividing it by windows, (b) Specific 
heat associated with the energy shown in the inset, which was obtained from the density in (a). 
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FIG. 3: CPDF for a = 1 and size L = 16, at three different temperatures. This figure clearly shows 
a first-order phase transition because of the double peaked form of the CPDF at the pseudo-critical 
temperature for this lattice size. 
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FIG. 4: Multi-range Wang-Landau results for the relative logarithm of the density of states, for 
fi = (\/eY, where i = 3. The blue curve represents the overlapped windows. 



11 




FIG. 5: (a) Specific heat for a = 0.5 and L = 16. The continuous line corresponds to WLS 
simulations, whereas red points to traditional Metropolis ones, (b) CPDF at the pseudo-critical 
temperature. 




FIG. 6: (a) Specific heat for a = 0.25 and L = 16. (b) CPDF at the corresponding pseudo-critical 
temperature. 
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FIG. 7: (a) Specific heat for a = 0.24 and L = 16. The black lines correspond to WLS results 
for different energy ranges, whereas red points to traditional Metropolis results, (b) CPDF at the 
corresponding pseudo-critical temperature for the specific heat peak obtained by WLS. The single 
peaked structure of the CPDF shows a second-order phase transition for the present size. 
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FIG. 9: Phase diagram of the model defined in Eq.(l) for L = 16. Circles represent estimated 
second-order critical points. Stars represent first-order critical points. The location of the critical 
end point is uncertain, so it must be close to the end of the arrow. The star at T = and a = 0.25 
is exactly located. All lines are just guides to the visualization. 
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